Auxiliary functions
get_filtered_genes_index <- function(scaled_tpm, meta_dt, min_count = 100, min_total = 150) {
design_matrix_groups <-
model.matrix(
~ condition:batch,
meta_dt %>%
column_to_rownames("sample")
)
keep <- edgeR::filterByExpr(as.matrix(scaled_tpm[, -c(1, 2)]),
design_matrix_groups,
min.count = min_count,
min.total.count = min_total
)
return(keep)
}
get_deseq_data <- function(target_tissue, design_matrix, gene_metadata, meta_dt, data_dir, tx2knownGene) {
samples_dir <- file.path(
data_dir,
meta_dt %>%
pull(name)
)
cdna_dir <- file.path(samples_dir)
cdna_files <- file.path(cdna_dir, "abundance.h5")
names(cdna_files) <-
meta_dt %>%
pull(sample)
# bias corrected counts without an offset. Corrects for potential differential isoform usage
txi_scaledTPM <- tximport(cdna_files,
type = "kallisto",
tx2gene = tx2knownGene,
countsFromAbundance = c("scaledTPM")
)
scaled_tpm <-
txi_scaledTPM$counts %>%
as_tibble(rownames = "gene_id")
scaled_tpm <-
left_join(scaled_tpm,
gene_metadata %>%
dplyr::select(gene_name, gene_id),
by = "gene_id"
) %>%
relocate(gene_name, .after = gene_id)
write_tsv(scaled_tpm, paste0(target_tissue, "_engraftment_scaled_tpm.tsv"))
# for deseq2
txi <- tximport(cdna_files,
type = "kallisto",
tx2gene = tx2knownGene,
countsFromAbundance = c("no")
)
dds <-
DESeqDataSetFromTximport(
txi,
meta_dt %>%
column_to_rownames("sample"),
design = design_matrix
)
keep <- get_filtered_genes_index(scaled_tpm, meta_dt)
table(keep)
dds <- dds[keep, ]
design_matrix_reduced <-
model.matrix(
~1,
meta_dt %>%
column_to_rownames("sample")
)
colnames(design_matrix_reduced) <- c("intercept")
dds <- DESeq(dds, test = "LRT", reduced = design_matrix_reduced)
counts_normalized <-
counts(dds, normalized = T) %>%
as_tibble(rownames = "gene_id") %>%
left_join(.,
gene_metadata %>%
dplyr::select(gene_name, gene_id),
by = "gene_id"
) %>%
relocate(gene_name, .after = gene_id)
write_tsv(counts_normalized, paste0(target_tissue, "_engraftment_counts_filtered_normalized.tsv"))
return(dds)
}
# We use apeglm shrinkage to preserve the size of large LFC and compute s-values (the estimated rate of false sign). Shrinkage is useful for ranking and visualization, without the need for arbitrary filters on low count genes. https://doi.org/10.1093/bioinformatics/bty895
# The local false sign rate FSR is defined as the posterior probability that the posterior mode (MAP) of beta (log2FC) is of a biologically significant effect size (i.e. abs(beta) > lfcThreshold )
get_lfc <- function(dds, coef_name, lfc_threshold, svalue_threshold) {
lfc <-
lfcShrink(dds,
coef = coef_name,
type = "apeglm",
svalue = T,
lfcThreshold = lfc_threshold
)
plotMA(lfc, ylim = c(-22, 8))
abline(h = c(-1, 1), col = "dodgerblue", lwd = 2)
lfc <-
lfc %>%
as_tibble() %>%
dplyr::select(-baseMean) %>%
mutate(
keep = ((abs(log2FoldChange) > lfc_threshold) & (svalue < svalue_threshold))
) %>%
set_names(c(
paste0(
c("lfc_", "lfc_se_", "svalue_", "keep_"),
coef_name
)
))
return(lfc)
}
get_lfc_estimates <- function(dds, gene_metadata, coef_names, lfc_threshold, svalue_threshold) {
lfcs <-
map_dfc(coef_names,
get_lfc,
dds = dds,
lfc_threshold = lfc_threshold,
svalue_threshold = svalue_threshold
)
#coef_names_batch <- c(paste0(c("intercept", coef_names), "_batch"), paste0("SE_", c("intercept", coef_names), "_batch"))
lfcs <-
bind_cols(
mcols(dds)[, c("baseMean", "baseVar", "dispOutlier", "maxCooks", "intercept", "SE_intercept", coef_names)] %>%
as_tibble(rownames = "gene_id") %>%
dplyr::select(gene_id, everything()),
lfcs
)
lfcs <-
lfcs %>%
dplyr::select(
gene_id, baseMean,
intercept, paste0("lfc_", coef_names),
paste0("svalue_", coef_names), everything()
)
lfcs <-
lfcs %>%
left_join(., gene_metadata, by = "gene_id") %>%
dplyr::select(gene_name, everything()) %>%
mutate(
gene_name =
scater::uniquifyFeatureNames(
ID = gene_id,
names = gene_name
)
)
return(lfcs)
}
plot_volcano <- function(data, xvar, yvar, title, lfc_threshold, svalue_thresh, figures_dir) {
gv <-
EnhancedVolcano(data,
lab = data %>% pull(gene_name),
xlab = bquote(~ italic(Moderated) ~ Log[2] ~ FC),
ylab = bquote(~ -Log[10] ~ italic(svalue)),
x = xvar,
y = yvar,
col = c("grey30", "forestgreen", "red2", "royalblue"),
pCutoff = svalue_thresh,
FCcutoff = lfc_threshold,
# legend = c("NS", "Log2FC", "svalue", "svalue & Log2FC"),
legendLabels = c(
"NS",
expression(Log[2] ~ FC),
"svalue",
expression(s - value ~ and ~ log[2] ~ FC)
)
)
ggsave(file.path(figures_dir, paste0(title, "_volcano.pdf")), gv)
return(gv)
}
plot_heatmap <- function(target_genes, title, meta_col, data, show_genes = F, figures_dir) {
df <-
left_join(target_genes %>%
dplyr::select(gene_id),
data,
by = "gene_id"
) %>%
dplyr::select(-gene_id) %>%
column_to_rownames("gene_name") %>%
as.matrix()
gh <-
pheatmap(df,
color = colorRampPalette(rev(RColorBrewer::brewer.pal(n = 11, name = "Spectral")))(100),
scale="row",
cluster_rows = T,
show_rownames = show_genes,
treeheight_row = 0,
fontsize_col = 4,
fontsize_row = 6,
angle_col = "45",
cluster_cols = FALSE,
annotation_col = meta_col
)
ggsave(file.path(figures_dir, paste0(title, "_heatmap_expression.pdf")), gh)
return(gh)
}
plot_pca <- function(vsd, title, figures_dir) {
pcaData <-
plotPCA(vsd,
intgroup = c("time", "age", "tissue", "batch", "donor"),
returnData = TRUE
)
percentVar <- round(100 * attr(pcaData, "percentVar"))
gg_plot <-
ggplot(
pcaData %>%
mutate(condition = paste(age, time, sep = "_")),
aes(x = PC1, y = PC2, color = time, shape = age, label = name)
) +
geom_point(size = 4) +
geom_line(aes(group = donor), colour = "grey") +
# coord_fixed() +
geom_text(check_overlap = T, angle = 0, size = 3, vjust = 3, nudge_y = 1.5) +
xlab(paste0("PC1: ", percentVar[1], "% variance")) +
ylab(paste0("PC2: ", percentVar[2], "% variance")) +
coord_fixed() +
ggtitle("PCA with VST data")
ggsave(file.path(figures_dir, paste0(title, "_pca_bulk_rnaseq.pdf")), gg_plot)
gg_plot
}
plot_gene <- function(gene, data, figures_dir) {
df_plot <-
data %>%
dplyr::filter(gene_id == gene)
gene_name <- df_plot$gene_name[1]
p <- ggplot(
df_plot,
aes(x = time, y = lognorm_counts)
) +
geom_point(aes(color = donor)) +
geom_line(aes(group = donor)) +
labs(title = gene_name) +
facet_grid(batch ~ age)
ggsave(file.path(figures_dir, paste0(gene_name, "_expression_lognorm_counts.pdf")), p)
return(p)
}
Get tx id versions
library(AnnotationHub)
Loading required package: BiocFileCache
Loading required package: dbplyr
Attaching package: 'dbplyr'
The following objects are masked from 'package:dplyr':
ident, sql
Attaching package: 'AnnotationHub'
The following object is masked from 'package:Biobase':
cache
ah <- AnnotationHub()
DEPRECATION: As of AnnotationHub (>2.23.2), default caching location has changed.
Problematic cache: /home/rick/.cache/AnnotationHub
See https://bioconductor.org/packages/devel/bioc/vignettes/AnnotationHub/inst/doc/TroubleshootingTheCache.html#default-caching-location-update
snapshotDate(): 2021-05-18
edb <- ah[["AH73905"]]
loading from cache
seqlevelsStyle(edb) <- "UCSC"
# standard_chroms <- standardChromosomes(edb, species = organism(edb))
tx_meta <-
transcripts(edb,
columns = c(
"tx_id_version",
"gene_id",
"gene_name",
"gene_id_version"
),
return.type = "DataFrame"
)
tx2knownGene <-
tx_meta %>%
as_tibble() %>%
dplyr::select(tx_id_version, gene_id)
meta_dt <-
read_tsv(file.path(metadata_filepath)) %>%
# arrange(tissue, age, time) %>%
mutate(
condition = paste(age, time, sep = "_"),
age = factor(age, levels = c("young", "aged")),
time = factor(time, levels = c("T0", "T21")),
batch = factor(batch, levels = c(1, 2))
)
Rows: 23 Columns: 8── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (7): name, sample, time, tissue, age, donor, type
dbl (1): batch
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
meta_dt
gene_metadata <-
read_csv(gene_metadata_filepath) %>%
dplyr::select(-gene_id_version)
Rows: 55487 Columns: 8── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (6): gene_name, gene_id, seqnames, gene_biotype, gene_id_version, description
dbl (2): gene_len, perc_gene_gc
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
gene_metadata <-
left_join(tx_meta %>% as_tibble() %>%
group_by(gene_id) %>%
dplyr::count(gene_name) %>%
dplyr::select(-n),
gene_metadata %>%
dplyr::select(-gene_name),
by = "gene_id"
) %>%
mutate(
gene_name =
scater::uniquifyFeatureNames(
ID = gene_id,
names = gene_name
)
)
gene_metadata
design_matrix <-
model.matrix(
~ time * age * batch,
meta_dt %>%
column_to_rownames("sample")
)
colnames(design_matrix) <- c(
"intercept",
"engraf",
"age",
"intercept_batch",
"niche",
"engraf_batch",
"age_batch",
"niche_batch"
)
design_matrix <- design_matrix[, c(1, 2, 3, 5, 4, 6, 7, 8)]
design_matrix
intercept engraf age niche intercept_batch engraf_batch age_batch niche_batch
yng_t0_01 1 0 0 0 0 0 0 0
yng_t0_02 1 0 0 0 0 0 0 0
yng_t0_03 1 0 0 0 0 0 0 0
yng_t0_04 1 0 0 0 1 0 0 0
yng_t0_05 1 0 0 0 1 0 0 0
yng_t0_06 1 0 0 0 1 0 0 0
yng_t21_01 1 1 0 0 0 0 0 0
yng_t21_02 1 1 0 0 0 0 0 0
yng_t21_04 1 1 0 0 1 1 0 0
yng_t21_05 1 1 0 0 1 1 0 0
yng_t21_06 1 1 0 0 1 1 0 0
aged_t0_01 1 0 1 0 0 0 0 0
aged_t0_02 1 0 1 0 0 0 0 0
aged_t0_03 1 0 1 0 1 0 1 0
aged_t0_04 1 0 1 0 1 0 1 0
aged_t0_05 1 0 1 0 1 0 1 0
aged_t0_06 1 0 1 0 1 0 1 0
aged_t21_01 1 1 1 1 0 0 0 0
aged_t21_02 1 1 1 1 0 0 0 0
aged_t21_03 1 1 1 1 1 1 1 1
aged_t21_04 1 1 1 1 1 1 1 1
aged_t21_05 1 1 1 1 1 1 1 1
aged_t21_06 1 1 1 1 1 1 1 1
Load data
dds <- get_deseq_data("leg", design_matrix, gene_metadata, meta_dt, data_dir, tx2knownGene)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
summarizing abundance
summarizing counts
summarizing length
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
summarizing abundance
summarizing counts
summarizing length
using counts and average transcript lengths from tximport
using supplied model matrix
estimating size factors
using 'avgTxLength' from assays(dds), correcting for library size
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
2 rows did not converge in beta, labelled in mcols(object)$fullBetaConv. Use larger maxit argument with nbinomLRT
dds
class: DESeqDataSet
dim: 17153 23
metadata(1): version
assays(6): counts avgTxLength ... H cooks
rownames(17153): ENSMUSG00000000001 ENSMUSG00000000028 ... ENSMUSG00000118487 ENSMUSG00000118488
rowData names(33): baseMean baseVar ... deviance maxCooks
colnames(23): yng_t0_01 yng_t0_02 ... aged_t21_05 aged_t21_06
colData names(8): name time ... batch condition
plotDispEsts(dds)

res <- results(dds)
W <- res$stat
maxCooks <- apply(assays(dds)[["cooks"]], 1, max)
idx <- !is.na(W)
plot(rank(W[idx]), maxCooks[idx],
xlab = "rank of Wald statistic",
ylab = "maximum Cook's distance per gene",
ylim = c(0, 5), cex = .4, col = rgb(0, 0, 0, .3)
)
m <- ncol(dds)
p <- 3
abline(h = qf(.99, p, m - p))

plot(res$baseMean + 1, -log10(res$pvalue),
log = "x", xlab = "mean of normalized counts",
ylab = expression(-log[10](pvalue)),
ylim = c(0, 30),
cex = .4, col = rgb(0, 0, 0, .3)
)

lfc_threshold <- 1
svalue_thresh <- 0.05
coef_names <- colnames(design_matrix)[-1]
lfcs <-
get_lfc_estimates(
dds,
gene_metadata,
coef_names,
lfc_threshold,
svalue_thresh
)
using 'apeglm' for LFC shrinkage. If used in published research, please cite:
Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
sequence count data: removing the noise and preserving large differences.
Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
sequence count data: removing the noise and preserving large differences.
Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
sequence count data: removing the noise and preserving large differences.
Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
sequence count data: removing the noise and preserving large differences.
Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
sequence count data: removing the noise and preserving large differences.
Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
sequence count data: removing the noise and preserving large differences.
Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)


lfcs
ggplot(lfcs, aes(x = lfc_niche_batch, y = lfc_age_batch)) +
geom_point(size = 1)

lfcs <-
lfcs %>%
mutate(
total_batch = abs(intercept_batch) + abs(niche_batch) + abs(age_batch) + abs(engraf_batch),
keep_niche_b = (keep_niche & !keep_niche_batch),
keep_age_b = (keep_age & !keep_age_batch),
keep_engraf_b = (keep_engraf & !keep_engraf_batch),
EAN = paste0(
as.integer(keep_engraf_b),
as.integer(keep_age_b),
as.integer(keep_niche_b)
),
)
table(lfcs$EAN)
000 001 010 011 100 101 110 111
15374 608 314 118 351 78 240 70
write_tsv(lfcs, "leg_two_batches_moderated_lfc_1_svalue_05_marked.txt")
summary_dt <-
colSums(assay(dds)) %>%
tibble::enframe(name = "sample", value = "counts")
summary_dt
scaled tpm
scaled_tpm <- read_tsv("leg_engraftment_scaled_tpm.tsv")
Rows: 36558 Columns: 25── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (2): gene_id, gene_name
dbl (23): yng_t0_01, yng_t0_02, yng_t0_03, yng_t0_04, yng_t0_05, yng_t0_06, yng_t21_01, yng_t21_02, yng_t21_04, yng_t21_05, yng_t21_06, aged_t0_01, aged_t0_02, aged_t0_03, aged_t0_04, aged_t0_05, aged_t0_06, aged_t21_01, aged_t21_02, aged_t2...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
counts_normalized <- read_tsv("leg_engraftment_counts_filtered_normalized.tsv")
Rows: 17153 Columns: 25── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (2): gene_id, gene_name
dbl (23): yng_t0_01, yng_t0_02, yng_t0_03, yng_t0_04, yng_t0_05, yng_t0_06, yng_t21_01, yng_t21_02, yng_t21_04, yng_t21_05, yng_t21_06, aged_t0_01, aged_t0_02, aged_t0_03, aged_t0_04, aged_t0_05, aged_t0_06, aged_t21_01, aged_t21_02, aged_t2...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
counts_normalized
counts_normalized_long <-
counts_normalized %>%
pivot_longer(cols = yng_t0_01:aged_t21_06, names_to = "sample", values_to = "counts") %>%
left_join(., meta_dt %>% dplyr::select(sample, time, age, donor, batch, condition), by = "sample") %>%
mutate(lognorm_counts = log1p(counts))
lfcs_hits <-
lfcs %>%
dplyr::filter(lfcs$EAN != "000")
lfcs_hits
targets <- c("Ccnd1", "Sparc", "Col3a1", "Mt2", "Ccl11", "Col4a2")
lfcs_genes <-
lfcs %>%
dplyr::filter(gene_name %in% targets)
lfcs_genes
lfcs_genes_dt<-
lfcs_genes %>%
dplyr::select(gene_name, EAN, intercept:lfc_niche, intercept_batch:niche_batch, svalue_engraf:svalue_niche, keep_engraf, keep_age, keep_niche, total_batch)
lfcs_genes_dt
figures_dir <- "./figures/hits"
map(lfcs_genes %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)
Saving 7 x 7 in image
[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]






Hits by group
target_genes_leg_001 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("001")) %>%
arrange(svalue_niche)
target_genes_leg_001
figures_dir <- "./figures/hits/leg_001"
pp1 <- map(target_genes_leg_001 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)
Saving 7 x 7 in image
target_genes_leg_011 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("011")) %>%
arrange(-baseMean)
target_genes_leg_011
figures_dir <- "./figures/hits/leg_011"
pp2 <- map(target_genes_leg_011 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)
Saving 7 x 7 in image
target_genes_leg_101 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("101"))
target_genes_leg_101
figures_dir <- "./figures/hits/leg_101"
pp3 <- map(target_genes_leg_101 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)
Saving 7 x 7 in image
target_genes_leg_111 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("111"))
target_genes_leg_111
figures_dir <- "./figures/hits/leg_111"
pp4 <- map(target_genes_leg_111 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)
Saving 7 x 7 in image
target_genes_leg_100 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("100"))
target_genes_leg_100
figures_dir <- "./figures/hits/leg_100"
pp5 <- map(target_genes_leg_100 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)
Saving 7 x 7 in image
target_genes_leg_010 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("010"))
target_genes_leg_010
figures_dir <- "./figures/hits/leg_010"
pp6 <- map(target_genes_leg_010 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)
Saving 7 x 7 in image
target_genes_leg_110 <-
lfcs_hits %>%
dplyr::filter(EAN %in% c("110"))
target_genes_leg_110
figures_dir <- "./figures/hits/leg_110"
pp7 <- map(target_genes_leg_110 %>% pull(gene_id), plot_gene, data = counts_normalized_long, figures_dir)
Saving 7 x 7 in image
svalue_thresh <- 0.05
figures_dir <- "./figures/volcanos"
plot_volcano(lfcs, "lfc_niche", "svalue_niche", "niche", lfc_threshold, svalue_thresh, figures_dir)
One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

plot_volcano(lfcs, "lfc_age", "svalue_age", "age", lfc_threshold, svalue_thresh, figures_dir)
Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

plot_volcano(lfcs, "lfc_engraf", "svalue_engraf", "engraf", lfc_threshold, svalue_thresh, figures_dir)
Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

svalue_thresh <- 0.05
figures_dir <- "./figures/volcanos"
plot_volcano(lfcs, "lfc_niche_batch", "svalue_niche_batch", "niche_batch", lfc_threshold, svalue_thresh, figures_dir)
Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

svalue_thresh <- 0.05
figures_dir <- "./figures/volcanos"
plot_volcano(lfcs, "lfc_age_batch", "svalue_age_batch", "age_batch", lfc_threshold, svalue_thresh, figures_dir)
Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

svalue_thresh <- 0.05
figures_dir <- "./figures/volcanos"
plot_volcano(lfcs, "lfc_engraf_batch", "svalue_engraf_batch", "engraf_batch", lfc_threshold, svalue_thresh, figures_dir)
Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

svalue_thresh <- 0.05
figures_dir <- "./figures/volcanos"
plot_volcano(lfcs, "lfc_intercept_batch", "svalue_intercept_batch", "intercept_batch", lfc_threshold, svalue_thresh, figures_dir)
One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...Ignoring unknown parameters: xlim, ylimSaving 7 x 7 in image

sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] parallel stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] AnnotationHub_3.0.2 BiocFileCache_2.0.0 dbplyr_2.1.1 edgeR_3.34.1 limma_3.48.3 RColorBrewer_1.1-2 PoiClaClu_1.0.2.1 scater_1.20.1
[9] scuttle_1.2.1 SingleCellExperiment_1.14.1 ensembldb_2.16.4 AnnotationFilter_1.16.0 GenomicFeatures_1.44.2 AnnotationDbi_1.54.1 pheatmap_1.0.12 EnhancedVolcano_1.10.0
[17] ggrepel_0.9.1 rhdf5_2.36.0 apeglm_1.14.0 glmpca_0.2.0 cowplot_1.1.1 vsn_3.60.0 DESeq2_1.32.0 SummarizedExperiment_1.22.0
[25] Biobase_2.52.0 MatrixGenerics_1.4.3 matrixStats_0.61.0 GenomicRanges_1.44.0 GenomeInfoDb_1.28.4 IRanges_2.26.0 S4Vectors_0.30.2 BiocGenerics_0.38.0
[33] tximport_1.20.0 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4 readr_2.1.2 tidyr_1.1.4 tibble_3.1.6
[41] ggplot2_3.3.5 tidyverse_1.3.1
loaded via a namespace (and not attached):
[1] utf8_1.2.2 tidyselect_1.1.1 RSQLite_2.2.9 grid_4.1.2 BiocParallel_1.26.2 munsell_0.5.0 ScaledMatrix_1.0.0
[8] ragg_1.2.1 preprocessCore_1.54.0 withr_2.4.3 colorspace_2.0-2 filelock_1.0.2 ggalt_0.4.0 knitr_1.37
[15] rstudioapi_0.13 Rttf2pt1_1.3.9 labeling_0.4.2 bbmle_1.0.24 GenomeInfoDbData_1.2.6 farver_2.1.0 bit64_4.0.5
[22] coda_0.19-4 vctrs_0.3.8 generics_0.1.1 xfun_0.29 R6_2.5.1 ggbeeswarm_0.6.0 rsvd_1.0.5
[29] locfit_1.5-9.4 bitops_1.0-7 rhdf5filters_1.4.0 cachem_1.0.6 DelayedArray_0.18.0 assertthat_0.2.1 vroom_1.5.7
[36] promises_1.2.0.1 BiocIO_1.2.0 scales_1.1.1 beeswarm_0.4.0 gtable_0.3.0 beachmat_2.8.1 ash_1.0-15
[43] affy_1.70.0 rlang_1.0.0 genefilter_1.74.1 systemfonts_1.0.3 splines_4.1.2 lazyeval_0.2.2 rtracklayer_1.52.1
[50] extrafontdb_1.0 broom_0.7.12 BiocManager_1.30.16 yaml_2.2.2 modelr_0.1.8 backports_1.4.1 httpuv_1.6.5
[57] extrafont_0.17 tools_4.1.2 affyio_1.62.0 ellipsis_0.3.2 jquerylib_0.1.4 Rcpp_1.0.8 plyr_1.8.6
[64] sparseMatrixStats_1.4.2 progress_1.2.2 zlibbioc_1.38.0 RCurl_1.98-1.5 prettyunits_1.1.1 viridis_0.6.2 haven_2.4.3
[71] fs_1.5.2 magrittr_2.0.2 reprex_2.0.1 mvtnorm_1.1-3 ProtGenerics_1.24.0 mime_0.12 hms_1.1.1
[78] evaluate_0.14 xtable_1.8-4 XML_3.99-0.8 emdbook_1.3.12 readxl_1.3.1 gridExtra_2.3 compiler_4.1.2
[85] biomaRt_2.48.3 bdsmatrix_1.3-4 maps_3.4.0 KernSmooth_2.23-20 crayon_1.4.2 htmltools_0.5.2 later_1.3.0
[92] tzdb_0.2.0 geneplotter_1.70.0 lubridate_1.8.0 DBI_1.1.2 proj4_1.0-11 MASS_7.3-55 rappdirs_0.3.3
[99] Matrix_1.4-0 cli_3.1.1 pkgconfig_2.0.3 GenomicAlignments_1.28.0 numDeriv_2016.8-1.1 xml2_1.3.3 annotate_1.70.0
[106] vipor_0.4.5 bslib_0.3.1 XVector_0.32.0 rvest_1.0.2 digest_0.6.29 Biostrings_2.60.2 rmarkdown_2.11
[113] cellranger_1.1.0 DelayedMatrixStats_1.14.3 restfulr_0.0.13 curl_4.3.2 shiny_1.7.1 Rsamtools_2.8.0 rjson_0.2.21
[120] lifecycle_1.0.1 jsonlite_1.7.3 Rhdf5lib_1.14.2 BiocNeighbors_1.10.0 viridisLite_0.4.0 fansi_1.0.2 pillar_1.6.5
[127] lattice_0.20-45 ggrastr_1.0.1 KEGGREST_1.32.0 fastmap_1.1.0 httr_1.4.2 survival_3.2-13 interactiveDisplayBase_1.30.0
[134] glue_1.6.1 png_0.1-7 BiocVersion_3.13.1 bit_4.0.4 stringi_1.7.6 sass_0.4.0 blob_1.2.2
[141] textshaping_0.3.6 BiocSingular_1.8.1 memoise_2.0.1 irlba_2.3.5